Galton Board

WLJS Version adapted from [Galton Board
2025, Gabriel Lemieux](https://www.wolframcloud.com/obj/gabriel.lemieux/Published/GaltonBoard.nb)
(* Initialization *)
ball = <|
  r -> 0.1,         (* Radius of the ball *)
  n -> 100          (* Number of simultaneous balls *)
|>;

peg = <|
  r -> 0.25         (* Radius of the pegs *)
|>;

sim = <|
  \[CapitalDelta]t -> 0.05,    (* Higher Value need to correct collision positions *)
  \[Epsilon] -> 0.4,           (* 1.0 means all energy back in the ball *)
  g -> {0, -9.81},             (* You can tilt the board *)
  n -> 1000                    (* Total Number of Balls *)
|>;

(* Initialize Position and Momentum Matrices *)
randomBall := {RandomReal[{-0.05, 0.05}], RandomReal[{0, 1}] + 10};
ball[pos] = Table[randomBall, ball[n]];
ball[mom] = Table[0, ball[n], 2];
ball[bins] = {};

(* Peg position *)
peg[pos] = Flatten[Table[{ii + 0.5 jj, jj}, {ii, -15, 10}, {jj, -1, 8}], 1];

(* Display function *)
positions = ball[pos];
Graphics[{
  {LightGray, EdgeForm[Gray], Disk[{#1, #2}, peg[r]]} & @@@ peg[pos],
  {PointSize[0.02], Opacity[0.8], Blue, Point[positions // Offload]},
  AnimationFrameListener[positions // Offload, "Event"->"frameTick"]
},
  AspectRatio -> 0.5,
  PlotRange -> {10 {-1, 1}, 10 {0, 1}},
  PlotRangeClipping -> True,
  ImageSize -> Large
]


EventHandler["frameTick", Function[Null,
If[Length[ball[pos]] > 0,

  (* Step *)
  (* Scheme: Backward Euler *)

  (* Next Step Momentum *)
  ball[mom] += Threaded[sim[\[CapitalDelta]t] sim[g]];

  (* Next Position Prediction *)
  next = ball[pos] + sim[\[CapitalDelta]t] ball[mom];

  (* Distance: Do not parallelize, it's already in the algo. *)
  dist = DistanceMatrix[next, peg[pos]];

  (* Collisions between current and future step *)
  col = Parallelize @ MapIndexed[
    {v, i} |-> {
      SelectFirst[v, (# < ball[r] + peg[r]) & -> "Index"], i[[1]]
    } /. ({Missing[_], _} -> Nothing),
    dist
  ];

  (* Correct Momentum to include collisions *)
  If[Length[col] > 1,
    {
      up, val
    } = Transpose[
      Parallelize @ MapApply[
        {p, b} |-> {
          b,
          sim[\[Epsilon]] * Norm[ball[mom][[b]]] Normalize[next[[b]] - peg[pos][[p]]]
        },
        col
      ]
    ];
    momt = ball[mom];
    momt[[up]] = val;
    ball[mom] = momt;
  ];

  (* Update pos with corrected momentum *)
  ball[pos] += sim[\[CapitalDelta]t] ball[mom];

  (* Filter through the balls, remove/reset and register those exiting *)
  exited = Select[ball[pos], #[[2]] < 0 & -> "Index"];
  ingame = Select[ball[pos], #[[2]] > 0 & -> "Index"];
  If[Length[exited] > 0,
    ball[bins] = Join[ball[bins], ball[pos][[exited, 1]]]
  ];

  (* Remove them from the pos and mom matrices *) 
  ball[pos] = ball[pos][[ingame]];
  ball[mom] = ball[mom][[ingame]];

  (* Insert new ball *)
  left = sim[n] - Length[ball[bins]] - Length[ball[pos]];
  toAdd = Min[left, ball[n] - Length[ball[pos]]];
  Do[
    AppendTo[ball[pos], randomBall];
    AppendTo[ball[mom], {0, 0}],
    toAdd
  ];

  positions = ball[pos];
]]];